npj Systems Biology and Applications

Cellular responses to perturbations and drugs are determined by interconnected networks, rather than linear pathways. Individually, the JNK, p38 and p53 stress and DNA-damage response networks are well understood and regulate critical cell-fate decisions, such as apoptosis, in response to many chemotherapeutical agents, such as doxorubicin. To better understand how interactions between these pathways determine the dynamic behaviour of the entire network, we constructed a data-driven mathematical model. This model contains mechanistic details about the kinase cascades that activate JNK, p38, AKT and p53, and free parameters that describe possible interactions between these pathways. Fitting this model to experimental time-course perturbation data (five time-courses with six time-points under five different conditions), identified specific network interactions that can explain the observed network responses. JNK emerged as an important control node. JNK exhibited a positive feedback loop, was tightly controlled by negative feedback and crosstalk from p38 and AKT, respectively, and was the strongest activator of p53. Compared to static network reconstruction methods, such as modular response analysis, the model-based approach identifies biochemical mechanisms and explains the dynamic control of cell signalling.

Difficult-to-treat (D2T) rheumatic disease affects approximately 12% of rheumatoid arthritis patients and resists sequential biologic therapy, yet no mechanistic model explains this refractoriness as a system-level phenomenon. Here we present the 3-Axis Integrative Framework (3-AIF), a six-variable ordinary differential equation system integrating mucosal tolerance, energy-gated neuroimmune danger sensing, and integrated stress response pathways coupled through Hill-function metabolic gating. Stability analysis reveals bistable dynamics with two co-existing attractors separated by a saddle point. Bifurcation analysis demonstrates fold catastrophe with hysteresis: recovery requires greater therapeutic effort than disease prevention. Sensitivity analysis identifies three dominant parameters mapping to neuroimmune activation, energy drain, and recovery capacity. Cross-disease transcriptomic consistency analysis across six public datasets (n=310, five rheumatic diseases, four tissue types) reveals compartment-specific axis dysregulation -- circulating cells show integrated stress response activation while target tissues show pathway exhaustion -- and disease-specific axis dominance patterns consistent with model predictions.

Systems biology aims to develop predictive models that connect molecular mechanisms to cellular behavior. Genome-scale metabolic models are among the most widely used frameworks for integrating stoichiometric, thermodynamic, and omics-derived information to predict feasible metabolic phenotypes. However, cellular metabolism operates on timescales governed by enzyme kinetics and by the relationship between metabolic fluxes and metabolite pool sizes. In steady-state metabolic models, this relationship can be expressed in terms of metabolite turnover rates, defined as flux-to-pool-size ratios that quantify how rapidly metabolite pools are renewed. As a result, physiologically consistent steady-state solutions should not only satisfy mass-balance and thermodynamic constraints but also exhibit turnover rates consistent with enzyme-mediated cellular dynamics. Current constraint-based approaches can admit many steady-state flux-concentration states that do not account for turnover rates, resulting in phenotypes incompatible with realistic metabolic dynamics, even when multiple types of data are imposed. Here, we present METEOR-K, an optimization framework that links steady-state metabolic fluxes to metabolite concentrations via turnover rate constraints to identify dynamically plausible flux-concentration reference states. Because these constraints reshape the feasible solution space, we also introduce turnover-rate-aware sampling strategies to efficiently explore the resulting feasible region. We applied METEOR-K to models of increasing scope and scale, including a reduced glycolysis pathway, anaerobic E. coli, and near-genome-scale ovarian cancer models. METEOR-K narrowed the admissible steady-state solution space, reduced uncertainty in feasible flux-concentration states, and improved local dynamic behavior. In nonlinear ODE simulations of bioreactor cultivation and drug-response scenarios, METEOR-K-derived states produced intracellular response times compatible with growth-supporting metabolic operation and perturbation recovery. Overall, these results establish metabolite turnover rates as scalable biophysical constraints that improve the physiological consistency of steady-state metabolic modeling. Because turnover rates encode flux-to-pool-size timescale constraints, METEOR-K moves part of physiological-consistency assessment upstream of kinetic parameterization, yielding better-suited flux-concentration reference states for kinetic modeling and dynamic prediction.

Treatment with drug combinations carries great promise for personalized therapy. We have previously shown that drug synergies targeting cancer can manually be identified based on a logical framework. We now demonstrate how automated adjustments of model topology and logic equations can greatly reduce the workload traditionally associated with logical model optimization. Our methodology allows the exploration of larger model ensembles that all obey a set of observations. We benchmark synergy predictions against a dataset of 153 targeted drug combinations. We show that well-performing manual models faithfully represent measured biomarker data and that their performance can be outmatched by automated parameterization using a genetic algorithm. The predictive performance of a curated model is strongly affected by simulated curation errors, while data-guided deletion of a small subset of edges can improve prediction quality. With correct topology we find some tolerance to simulated errors in the biomarker calibration data. With our framework we predict the synergy of joint inhibition of PI3K and TAK1, and further substantiate this prediction with observation in cancer cell cultures and in xenograft experiments.

Type 1 diabetes mellitus (T1D) is a systemic disease triggered by a local autoimmune inflammatory reaction in insulin-producing cells that disrupts the glucose-insulin-glucagon system and induces organ-wide, long-term effects on glycolytic and nonglycolytic processes. Mathematical modeling of the whole-body regulatory bihormonal system has helped to identify intervention points to ensure better control of T1D but was limited to a coarse-grained representation of metabolism. To extend the depiction of T1D, we developed a whole-body model using a novel integrative modeling framework that links organ-specific regulation and metabolism. The developed framework allowed the correct prediction of disrupted metabolic processes in T1D, highlighted pathophysiological processes common with neurodegenerative disorders, and suggested calcium channel blockers as potential adjuvants for diabetes control. Additionally, the model predicted the occurrence of insulin-dependent rewiring of interorgan crosstalk. Moreover, a simulation of a population of virtual patients allowed an assessment of the impact of inter and intraindividual variability on insulin treatment and the implications for clinical outcomes. In particular, GLUT4 was suggested as a potential pharmacogenomic regulator of intraindividual insulin efficacy. Taken together, the organ-resolved, dynamic model may pave the way for a better understanding of human pathology and model-based design of precise allopathic strategies.

The pancreatic islet (islet of Langerhans) is a mini-organ comprising several thousand endocrine cells, functioning jointly to maintain normoglycemia. Cellular networks within an islet were shown to influence its function in health and disease, but there are major gaps in our quantitative understanding of such architecture-function relations. Comprehensive modeling of an islet architecture and function requires the integration of vast amounts of information obtained through different experimental and theoretical approaches. To address this challenge, our lab has recently developed Bayesian metamodeling, a general approach for modeling complex systems by integrating heterogeneous input models. Here, we further developed metamodeling and applied it to construct a metamodel of a pancreatic islet. The metamodel relates islet architecture and function by combining a Monte-Carlo model of architecture trained on islet imaging data; and an ordinary differential equations (ODEs) mathematical model of function trained on calcium imaging, hormone imaging, and electrophysiological data. These input models are converted to a standardized statistical representation relying on Probabilistic Graphical Models; coupled by modeling their mutual relations with the physical world; and finally, harmonized through backpropagation. We validate the metamodel using existing data and use it to derive a testable hypothesis regarding the functional effect of varying intercellular connections. Since metamodeling currently requires substantial expert intervention, we also develop an automation tool for converting models to PGMs (step I) using feedforward neural networks. This automation is a first step towards automating the entire metamodeling process, working towards collaborative science through sharing of expertise, resources, data, and models.

The ability to reliably predict and infer cellular responses to environmental exposures would offer a major advance in the investigation of immune regulation in health and disease. One possible approach is the use of in silico modelling. Design of such a mathematical kinetic model would be based on existing knowledge of a biological system and utilise a partial data set to parameterise. However, the process of parameter estimation, key for the accuracy of the model, is difficult to conduct by hand, and thus a computational alternative is necessary. We report the utility of Genetic Algorithm with Rank Selection (GARS) as a parameter estimation tool on multiple biological models, including heat shock, signal transduction via ERK, circadian rhythm and NF{kappa}B systems, where it showed strong accuracy and superiority to the Extended Kalman Filter method, Algebraic Difference Equations, and MATLAB fminsearch approaches. GARS parameter estimation is a valuable tool for biological data because it reliably infers system behaviour from partial data sets, allowing for the prediction of cellular responses to environmental exposures.

Constraint-based metabolic modeling is a powerful way to study the mechanistic basis of cellular states and disease, but effective use demands substantial computational expertise and careful coordination of multi-step analyses. We developed MechAInistic to lower this barrier enabling researchers to ask complex biological questions in natural language. MechAInistic is a multi-agent system harnessing large language models organized around an Architect-Reviewer pattern that that converts a natural-language question into an executable, model-grounded workflow and produces a structured report. It supports pathway comparison, perturbation analysis, drug-target exploration, and literature interpretation across healthy and disease paired states. We evaluated MechAInistics therapeutic hypothesis generation using two immune-cell use-cases. For rheumatoid arthritis/healthy Naive B models, it identified mitochondrial metabolic rewiring and nominated Devimistat/CPI-613 as an investigational OGDH-centered hypothesis. In CD4+ Th17 multiple sclerosis/healthy models, the workflow identified NADP-dependent isocitrate dehydrogenase as the optimal target and proposed Ivosidenib as an FDA-approved repurposing candidate. GRAPHICAL ABSTRACT O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=83 SRC="FIGDIR/small/723319v1_ufig1.gif" ALT="Figure 1"> View larger version (19K): org.highwire.dtl.DTLVardef@1b5c1d1org.highwire.dtl.DTLVardef@1c798cforg.highwire.dtl.DTLVardef@10161d3org.highwire.dtl.DTLVardef@1bd7dce_HPS_FORMAT_FIGEXP M_FIG C_FIG

The immune system is a complex and dynamic network, crucial for combating infections and maintaining health. Developing a comprehensive digital twin of the immune system requires incorporating essential cellular components and their interactions. This study presents the first blueprint for an immune system digital twin, consisting of a comprehensive and simulatable mechanistic model. It integrates 51 innate and adaptive immune cells, 37 secretory factors, and 11 disease conditions, providing the foundation for developing a multi-scale model. The cellular-level model demonstrates its potential in characterizing immune responses to various single and combinatorial disease conditions. By making the model available in easy-to-use formats directly in the Cell Collective platform, the community can easily and further expand it. This blueprint represents a significant step towards developing general-purpose immune digital twins, with far-reaching implications for the future of digital twin technology in life sciences and healthcare, advancing patient care, and accelerating precision medicine.

The proper balance and transition between cellular quiescence and proliferation are critical to tissue homeostasis, repair, and regeneration. The likelihood of quiescence-to-proliferation transition is inversely correlated with quiescence depth, and deep quiescent cells are less likely to exit quiescence and reenter the cell cycle than shallow quiescent cells. The regulatory mechanisms of quiescence depth are poorly understood but essential for developing strategies against hypo- or hyper-proliferation diseases such as aging and cancer. Our earlier studies have demonstrated that the activation threshold of the bistable Rb-E2F gene network switch (ThE2F) controls quiescence depth. We have also identified coarse- and fine-tuning ThE2F modulators in rat embryonic fibroblasts. To examine whether other quiescent cells (including most adult stem and progenitor cells) under different environmental conditions use the same or different modulators of quiescence depth, here we studied the behaviors of 30,000 theoretical quiescent cell models that each support a functional Rb-E2F bistable switch with a unique parameter set. We found that although the vastly heterogeneous quiescent cell models exhibited no apparent parameter patterns, they converged at two alternative groups of strong quiescence-depth modulators (G1 cyclin/cdk-related and Rb/E2F complex-related). Our further machine learning (decision tree) analysis suggested that the Rb protein level and dephosphorylation rate in quiescent cells determine which modulator group to use to regulate quiescence depth.

Macrophages are crucial immune regulators as they can either trigger or resolve inflammation. These properties rely on defined inflammatory states and make macrophages valuable therapeutic targets. Identification of stable steady states and bistability in inflammatory transitions provides deeper insights into immune regulation and facilitates the development of novel therapeutic strategies. We present a multiomic and systems biology approach for the top-down identification of bistable circuits in human macrophages polarized towards pro- and anti-inflammatory phenotypes. Using differential gene expression profiles, we identified three criteria to suspect bistability: two potential attractors, a hysteresis behavior between their transitions, and the presence of potential modules of coregulated genes. This was further confirmed by proteomics data pointing to mutually exclusive and time-dependent profiles of gene expression. By network simplification and creation of a novel pipeline for parameter estimation in bistable models, we obtained a minimal model of inflammatory transitions in which we identified ultrasensitivity and hysteresis. Our minimal model genes establish a regulatory circuit switching miR-155 expression, which in turn regulates the expression of inflammatory marker genes during inflammatory transitions.

Symmetry refers to properties that remain invariant upon mathematical transformations. The principles of symmetry have guided numerous important discoveries in physics and chemistry but not in biology and medicine. Here, we aim to explore the presence of symmetry relationships at the gene expression level as a mean to distinguish between healthy and disease states. We deployed Learning-Based Invariant Feature Engineering - LIFE, a hybrid machine learning approach implemented with two symmetric invariant feature functions (IFFs) to identify Invariant Feature Genes (IFGs), which are gene pairs whose IFF single-value outputs remain invariant across individual samples in a given biological phenotype. Our multiclass classification results across the transcriptomes of 25 normal organs, 25 cancer types, and blood samples obtained from 4 different types of neurodegenerative diseases revealed the presence of unique phenotype-specific IFGs. We constructed networks using these IFGs (IF-Nets) and intriguingly, we demonstrated that the hubs could serve as information encoders, capable of reconstructing sample-wise expression values in relation to their counterpart genes. More importantly, we found that hubs of cancer IF-Nets were enriched with both approved and clinical trial drugs, highlighting "symmetry breaking" as a novel approach for treating diseases.

Yes-associated protein (YAP) and its homolog TAZ are transducers of several biochemical and biomechanical signals, serving to integrate multiplexed inputs from the microenvironment into higher-level cellular functions such as proliferation, differentiation, apoptosis, migration, and hemostasis. Emerging evidence suggests that Ca2+ is a key second messenger that closely connects microenvironmental input signals and YAP/TAZ regulation. However, studies that directly modulate Ca2+ have reported contradictory YAP/TAZ responses: In some studies, a reduction in Ca2+ influx increases the activity of YAP/TAZ, while in others, an increase in Ca2+ influx activates YAP/TAZ. Importantly, Ca2+ and YAP/TAZ exhibit distinct spatiotemporal dynamics, making it difficult to unravel their connections from a purely experimental approach. In this study, we developed a network model of Ca2+-mediated YAP/TAZ signaling to investigate how temporal dynamics and crosstalk of signaling pathways interacting with Ca2+ can alter YAP/TAZ response, as observed in experiments. By including six signaling modules (e.g., GPCR, IP3-Ca2+, Kinases, RhoA, F-actin, and Hippo-YAP/TAZ) that interact with Ca2+, we investigated both transient and steady-state cell response to Angiotensin II and thapsigargin stimuli. The model predicts stimuli, Ca2+ transient, and frequency-dependent relationships between Ca2+ and YAP/TAZ primarily mediated by signaling species like cPKC, DAG, CaMKII, and F-actin. Model results illustrate the role of Ca2+ dynamics and CaMKII bistable response in switching the direction of changes in Ca2+-induced YAP/TAZ activity for different stimuli. Frequency-dependent YAP/TAZ response revealed the competition between upstream regulators of LATS1/2, leading to the YAP/TAZ non-monotonic response to periodic GPCR stimulation. This study provides new insights into the underlying mechanisms responsible for the controversial Ca2+-YAP/TAZ relationship observed in experiments.

Modeling of systems for which data is limited often leads to underdetermined model identification problems, where multiple candidate models are equally adherent to data. In such situations additional optimality criteria are useful in model selection apart from the conventional minimization of error and model complexity. This work presents the attractor landscape as a domain for novel model selection criteria, where the number and location of attractors impact desirability. A set of candidate models describing immune response dynamics to SARS-CoV infection is used as an example for model selection based on features of the attractor landscape. Using this selection criteria, the initial set of 18 models is ranked and reduced to 7 models that have a composite objective value with a p-value < 0.05. Additionally, the impact of pharmacologically induced remolding of the attractor landscape is presented.

Understanding how distinct neuronal subtypes contribute to Alzheimers disease (AD) pathology remains a major challenge. Patient-derived induced pluripotent stem cell (iPSC) studies have shown neuronal subtype-specific molecular and pathological signatures, yet the underlying metabolic shifts driving this selective vulnerability are not completely understood. Here we present iNeuron-GEM, the first manually curated, genome-scale metabolic network of human neurons that integrates transcriptomic and metabolic knowledge to resolve subtype-specific metabolic states. By coupling iNeuron-GEM with single nucleus RNA sequencing data from post-mortem human cohort studies, ROSMAP and SEA-AD, we capture neuronal subtype-specific metabolic features and fluxes and identify perturbations in lipid and energy metabolism across excitatory and inhibitory neurons. Integrative analysis with NPS-AD data shows overlapping metabolic disruptions in AD and schizophrenia (SCZ), suggesting shared molecular vulnerabilities between neurodegenerative and neuropsychiatric disorders. We also developed a computational pipeline to infer transcriptional regulation of metabolic pathways and identify NR6A1 and NR3C1 as important regulators of lipid dysregulation in AD neurons. Our study establishes iNeuron-GEM as a framework to identify neuronal subtype-specific metabolic vulnerabilities in complex brain disorders.

Computational models of signaling networks provide frameworks for predicting how molecular cues guide cell decisions. But they are typically limited by manual curation from incomplete literature. Here, we test whether general-purpose large language models (LLMs) generate accurate models of signaling networks. We find that general purpose LLMs generate 24-58% of the reactions of literature-curated networks for cardiomyocyte hypertrophy, myofibroblast activation, and mechano-signaling, and predicting network responses to perturbations with accuracies of 5-26%. While current general-purpose LLMs generate signaling networks with limited accuracy, this study provides a pipeline and benchmarks to guide future improvements.

Although endocrine therapies and Cdk4/6 inhibitors have produced significantly improved outcomes for patients with estrogen receptor positive (ER+) breast cancer, continuous application of these drugs often results in resistance. We hypothesized that cancer cells acquiring drug resistance might increase their dependency on negative regulators of the cell cycle. Therefore, we investigated the effect of inhibiting WEE1 on delaying the development of resistance to palbociclib and fulvestrant. We treated ER+ MCF7 breast cancer cells with palbociclib alternating with a combination of fulvestrant and a WEE1 inhibitor AZD1775 for 12 months. We found that the alternating treatment prevented the development of drug resistance to palbociclib and fulvestrant compared to monotherapies. Furthermore, we developed a mathematical model that can simulate cell proliferation under monotherapy, combination or alternating drug treatments. Finally, we showed that the mathematical model can be used to minimize the number of fulvestrant plus AZD1775 treatment periods while maintaining its efficacy.

Biological networks that are formed through amalgamation of signaling pathways include recurrent configurations called network motifs. These statistically over-represented subgraphs are often formed through interconnected enzyme-substrate relationships that are known to result in highly dynamic downstream behavior, including oscillatory output. Such signals are abundant in biology: heartbeats, circadian rhythms, and cell cycles all exhibit characteristic frequencies. Though there has been great emphasis on how oscillations can be generated through network dynamics, little is known about oscillatory information processing or transduction capacity of network motifs. We employ ordinary differential equations-based dynamical modeling to understand how different network topologies impact oscillatory signal propagation through a multi-enzyme network. We model enzyme-substrate interactions of 20 commonly observed motifs using Michaelis-Menten kinetics. We then perform deterministic Monte Carlo simulations using a range of biologically relevant enzymatic parameters and input frequencies. From these simulations, we quantify signal propagation characteristics using cluster analysis, categorize different motif responses based on output characteristics, and identify potential mechanisms for oscillatory signal processing using parameter sensitivity analysis. We see that the input-output responses depend on network topology and enzyme kinetic parameters. Enzymatic motifs show median oscillatory suppression of 30-135 decibels, with three-node coherent feedforward loops showing the lowest propensity for oscillatory signal suppression. Motifs that contained negative feedback or four-node coherent feedforward loops had the biggest potential to act as AC-to-DC converters, translating oscillatory input signals into transient impulses or sustained continuous outputs, respectively. We conclude that enzyme networks can process and decode information within oscillatory inputs in a frequency- and network-dependent manner.

We propose a method for creating personalized regulatory networks, enabling the development of virtual avatars for cancer patients, with patient-derived xenograft (PDX) models as a test case. Starting from a Prior Knowledge Network (PKN) based on the hallmarks of cancer, we constructed gene networks that are contextualized to each sample by integrating sample-specific gene expression data. These networks were optimized using a genetic algorithm to align with individual molecular profiles, focusing on key cancer-related processes. Following network optimization, we employed Graph Convolutional Networks (GCNs) to classify samples based the structures and interactions of their individualized network models and molecular profiles. This personalized approach provides insights into drug responses and helps predict treatment outcomes, offering a path toward more targeted cancer therapies. Author summaryCancer treatment can be more effective when therapies are personalized to each patients unique molecular profile. In this study, we introduce a method to create virtual avatars of cancer patients by personalizing regulatory networks using patient-derived xenograft (PDX) models as a proof of concept. Starting from known cancer hallmarks, we developed individualized gene networks for each sample by leveraging their specific gene expression data. These networks were refined with an optimization process to match the distinct molecular characteristics of each sample. By applying advanced machine learning, specifically Graph Convolutional Networks (GCNs), we classified these personalized models to better understand likely drug responses and predict treatment outcomes. This approach brings us closer to tailoring cancer therapies to individual patients, potentially improving treatment success by targeting key cancer pathways unique to each person.

Protein tyrosine phosphatase receptor J (PTPRJ) restrains cell proliferation and migration by dephosphorylating receptor tyrosine kinases (RTKs) including the epidermal growth factor receptor (EGFR). PTPRJ is a purported tumor suppressor, and alterations to its expression and/or function are associated with colorectal, breast, lung, and other cancers. While there is interest in controlling PTPRJ-regulated phenotypes, efforts are limited by the complexity of PTPRJ-mediated signaling. PTPRJ dephosphorylates multiple RTKs, and the degree to which PTPRJ control of signaling and phenotypes depends on local cellular RTK activation profiles is unknown. To probe the context dependence of PTPRJ signaling regulation, we collected signaling measurements across 16 pathway nodes at two time points in a panel of HSC3 carcinoma cells engineered with different PTPRJ expression profiles. Cells were treated with three different RTK ligands, and paired phenotype measurements (viability, wound healing, xCELLigence cell index) were made. Partial least squares regression models were developed to predict relationships between PTPRJ-regulated signaling pathways and cell phenotypes. The model effectively separated contributions to variance arising from the PTPRJ expression background and growth factor context. In testing model predictions, we demonstrated that PTPRJ suppressed MET-induced cell cell proliferation via regulation of a HER3/AKT signaling axis that stabilized PTPRJ expression through an unanticipated feedback mechanism. We also found that PTPRJ regulated HSC3 cell migration via JNK signaling that was preferentially activated by MET. Our results identify new regulatory nodes through which PTPRJ influences cancer cell phenotypes and demonstrates that these processes preferentially occur in the context of distinct RTK activation states.